*Decomposition*
clear

use "Results\tab_B3\Scenario_3\joining_couples.dta"

*Drop the benchmark categorization and rename the new one to make it fit with code below

drop ambition am_male am_female

rename ambition_type_sep ambition
rename am_sep_male am_male
rename am_sep_female am_female

cap drop ambition_type_k_4_s

gen ambition_type_k_4_s=ambition

*Keep only necessary variables
drop aegte_nr faelle_nr civst educ_eika final_educ individual grad_region wage_start_mean_ambition* wage_growth_ambition* wage_start_mean_ambition_s* wage_growth_ambition_s* married cohab educ_level_field

*Preparation

{
*get 1980 preferences*

*couple matrix

tab am_male am_female if aar==1980 & koen=="1", matcell(am_base) /*It does not matter if we condition on the husband or wife*/

*single men
tab ambition_type_k_4_s if koen=="1" & relationship==0 & aar==1980, matcell(am_single_m_base)

*single women
tab ambition_type_k_4_s if koen=="2" & relationship==0 & aar==1980, matcell(temp)
mat am_single_f_base=temp'

*preference matrix

mat am_pi_base=J(4,4,.)
forvalues i=1(1)4{
	forvalues j=1(1)4{
		mat am_pi_base[`i',`j']=am_base[`i',`j']/((am_single_m_base[`i',1]*am_single_f_base[1,`j'])^(1/2))
	}
}

}

*get 1980 marginals*
mat am_mar_m_base=J(4,1,.)
forvalues i=1(1)4{
mat am_mar_m_base[`i',1]=am_single_m_base[`i',1]+am_base[`i',1]+am_base[`i',2]+am_base[`i',3]+am_base[`i',4]
}


mat am_mar_f_base=J(1,4,.)
forvalues i=1(1)4{
mat am_mar_f_base[1,`i']=am_single_f_base[1,`i']+am_base[1,`i']+am_base[2,`i']+am_base[3,`i']+am_base[4,`i']
}




*get later preferences*

forvalues y=2018(1)2018{
tab am_male am_female if aar==`y' & koen=="1", matcell(am_`y') /*It does not matter if we condition on the husband or wife*/


*single men
tab ambition_type_k_4_s if koen=="1" & relationship==0 & aar==`y', matcell(am_single_m_`y')


*single women
tab ambition_type_k_4_s if koen=="2" & relationship==0 & aar==`y', matcell(temp)
mat am_single_f_`y'=temp'


*preference matrix

mat am_pi_`y'=J(4,4,.)
forvalues i=1(1)4{
	forvalues j=1(1)4{
		mat am_pi_`y'[`i',`j']=am_`y'[`i',`j']/((am_single_m_`y'[`i',1]*am_single_f_`y'[1,`j'])^(1/2))
	}
}



*get later marginals*
mat am_mar_m_`y'=J(4,1,.)
forvalues i=1(1)4{
mat am_mar_m_`y'[`i',1]=am_single_m_`y'[`i',1]+am_`y'[`i',1]+am_`y'[`i',2]+am_`y'[`i',3]+am_`y'[`i',4]
}


mat am_mar_f_`y'=J(1,4,.)
forvalues i=1(1)4{
mat am_mar_f_`y'[1,`i']=am_single_f_`y'[1,`i']+am_`y'[1,`i']+am_`y'[2,`i']+am_`y'[3,`i']+am_`y'[4,`i']
}



}

**Inequality**

keep if aar==1980 | aar==2018

gen couple_id_with_s=couple_id
replace couple_id_with_s=pnr if relationship==0
sort couple_id_with_s aar koen

gen individual_earnings=erhvervsindk_13
replace individual_earnings=0 if individual_earnings<0

by couple_id_with_s aar: egen household_earnings=sum(individual_earnings)

**ADDITION**
replace household_earnings=household_earnings/1.5 if relationship==1
**

keep if (koen=="1" & relationship==1) | relationship==0 

*Keep only necessary variables
drop erhvervsindk_13 couple_id couple_id_with_s individual_earnings

**Data**
{

*Data with singles*
foreach y in 1980 2018{
qui ineqdeco household_earnings if aar==`y'

sca gini_`y'=r(gini)
sca p90p50_`y'=r(p90p50)
sca p50p10_`y'=r(p50)/r(p10)
}

frame copy default results
frame change results
gen temp=.
collapse (first) temp, by(aar)
drop temp

gen gini=.
foreach i in 1980 2018{
replace gini=gini_`i' if aar==`i' 
}

gen p90p50=.
foreach i in 1980 2018{
replace p90p50=p90p50_`i' if aar==`i' 
}

gen p50p10=.
foreach i in 1980 2018{
replace p50p10=p50p10_`i' if aar==`i' 
}

frame change default


}



*Fixed household composition*

{
*Solve for 2018 single numbers with 1980 preferences and 1980 marginals (both)*

*ambition types

program nlfaq

syntax varlist(min=1 max=1) [if], at(name)

tempname m1 m2 m3 m4 f1 f2 f3 f4

scalar `m1' = `at'[1,1]
scalar `m2' = `at'[1,2]
scalar `m3' = `at'[1,3]
scalar `m4' = `at'[1,4]
scalar `f1' = `at'[1,5]
scalar `f2' = `at'[1,6]
scalar `f3' = `at'[1,7]
scalar `f4' = `at'[1,8]

tempvar yh
gen double `yh'= am_pi_base[1,1]*((`m1'*`f1')^(1/2))+am_pi_base[1,2]*((`m1'*`f2')^(1/2))+am_pi_base[1,3]*((`m1'*`f3')^(1/2))+am_pi_base[1,4]*((`m1'*`f4')^(1/2))+`m1'-am_mar_m_base[1,1]+1 in 1
replace `yh'=am_pi_base[2,1]*((`m2'*`f1')^(1/2))+am_pi_base[2,2]*((`m2'*`f2')^(1/2))+am_pi_base[2,3]*((`m2'*`f3')^(1/2))+am_pi_base[2,4]*((`m2'*`f4')^(1/2))+`m2'-am_mar_m_base[2,1]  in 2
replace `yh'=am_pi_base[3,1]*((`m3'*`f1')^(1/2))+am_pi_base[3,2]*((`m3'*`f2')^(1/2))+am_pi_base[3,3]*((`m3'*`f3')^(1/2))+am_pi_base[3,4]*((`m3'*`f4')^(1/2))+`m3'-am_mar_m_base[3,1] in 3
replace `yh'=am_pi_base[4,1]*((`m4'*`f1')^(1/2))+am_pi_base[4,2]*((`m4'*`f2')^(1/2))+am_pi_base[4,3]*((`m4'*`f3')^(1/2))+am_pi_base[4,4]*((`m4'*`f4')^(1/2))+`m4'-am_mar_m_base[4,1] in 4
replace `yh'=am_pi_base[1,1]*((`m1'*`f1')^(1/2))+am_pi_base[2,1]*((`m2'*`f1')^(1/2))+am_pi_base[3,1]*((`m3'*`f1')^(1/2))+am_pi_base[4,1]*((`m4'*`f1')^(1/2))+`f1'-am_mar_f_base[1,1] in 5
replace `yh'=am_pi_base[1,2]*((`m1'*`f2')^(1/2))+am_pi_base[2,2]*((`m2'*`f2')^(1/2))+am_pi_base[3,2]*((`m3'*`f2')^(1/2))+am_pi_base[4,2]*((`m4'*`f2')^(1/2))+`f2'-am_mar_f_base[1,2] in 6
replace `yh'=am_pi_base[1,3]*((`m1'*`f3')^(1/2))+am_pi_base[2,3]*((`m2'*`f3')^(1/2))+am_pi_base[3,3]*((`m3'*`f3')^(1/2))+am_pi_base[4,3]*((`m4'*`f3')^(1/2))+`f3'-am_mar_f_base[1,3] in 7
replace `yh'=am_pi_base[1,4]*((`m1'*`f4')^(1/2))+am_pi_base[2,4]*((`m2'*`f4')^(1/2))+am_pi_base[3,4]*((`m3'*`f4')^(1/2))+am_pi_base[4,4]*((`m4'*`f4')^(1/2))+`f4'-am_mar_f_base[1,4] in 8


replace `varlist'=`yh'
end


frame create temp
frame change temp

clear
set obs 8

gen y = 0
replace y=1 in 1

nl faq @ y, parameters(m1 m2 m3 m4 f1 f2 f3 f4) initial(m1 1 m2 1 m3 1 m4 1 f1 1 f2 1 f3 1 f4 1)

*get counter factual singles
mat am_single_m_cf_sort_2018=J(4,1,.)
mat am_single_m_cf_sort_2018[1,1]=[m1]_b[_cons]
mat am_single_m_cf_sort_2018[2,1]=[m2]_b[_cons]
mat am_single_m_cf_sort_2018[3,1]=[m3]_b[_cons]
mat am_single_m_cf_sort_2018[4,1]=[m4]_b[_cons]

mat am_single_f_cf_sort_2018=J(1,4,.)
mat am_single_f_cf_sort_2018[1,1]=[f1]_b[_cons]
mat am_single_f_cf_sort_2018[1,2]=[f2]_b[_cons]
mat am_single_f_cf_sort_2018[1,3]=[f3]_b[_cons]
mat am_single_f_cf_sort_2018[1,4]=[f4]_b[_cons]


frame change default
frame drop temp
prog drop nlfaq



*get counter factual couple numbers*

mat am_cf_sort_2018=J(4,4,.)

forvalues i=1(1)4{
	forvalues j=1(1)4{
		mat am_cf_sort_2018[`i',`j']=am_pi_base[`i',`j']*((am_single_m_cf_sort_2018[`i',1]*am_single_f_cf_sort_2018[1,`j'])^(1/2))
	}
}

}


*obtain weights - ambition

frame change results

gen gini_cf_sort=.
gen p90p50_cf_sort=.
gen p50p10_cf_sort=.

gen gini_m_cf_sort=.
gen p90p50_m_cf_sort=.
gen p50p10_m_cf_sort=.

frame change default

forvalues y=2018(1)2018{
	
	
*numerator 1
mat am_num_1_couples=J(4,4,.) /*The prob. of counterfactual given couple cell x,x*/

forvalues i=1(1)4{
	forvalues j=1(1)4{
		mat am_num_1_couples[`i',`j']=am_cf_sort_2018[`i',`j']/(am_cf_sort_2018[`i',`j']+am_`y'[`i',`j'])
	}
}

mat am_num_1_single_m=J(4,1,.) /*The prob. of counterfactual given single men cell x*/

forvalues i=1(1)4{
	mat am_num_1_single_m[`i',1]=am_single_m_cf_sort_2018[`i',1]/(am_single_m_cf_sort_2018[`i',1]+am_single_m_`y'[`i',1])
}

mat am_num_1_single_f=J(1,4,.) /*The prob. of counterfactual given single women cell x*/

forvalues i=1(1)4{
	mat am_num_1_single_f[1,`i']=am_single_f_cf_sort_2018[1,`i']/(am_single_f_cf_sort_2018[1,`i']+am_single_f_`y'[1,`i'])
}

*denominator 1
mat am_de_1_couples=J(4,4,.) /*The prob. of factual given couple cell x,x*/

forvalues i=1(1)4{
	forvalues j=1(1)4{
		mat am_de_1_couples[`i',`j']=am_`y'[`i',`j']/(am_cf_sort_2018[`i',`j']+am_`y'[`i',`j'])
	}
}

mat am_de_1_single_m=J(4,1,.) /*The prob. of factual given single men cell x*/

forvalues i=1(1)4{
	mat am_de_1_single_m[`i',1]=am_single_m_`y'[`i',1]/(am_single_m_cf_sort_2018[`i',1]+am_single_m_`y'[`i',1])
}

mat am_de_1_single_f=J(1,4,.) /*The prob. of counterfactual given single women cell x*/

forvalues i=1(1)4{
	mat am_de_1_single_f[1,`i']=am_single_f_`y'[1,`i']/(am_single_f_cf_sort_2018[1,`i']+am_single_f_`y'[1,`i'])
}

*numerator 2 /*The prob. of being in factual sample*/

sca number_factual=am_`y'[1,1]+am_`y'[1,2]+am_`y'[1,3]+am_`y'[1,4] ///
	+am_`y'[2,1]+am_`y'[2,2]+am_`y'[2,3]+am_`y'[2,4] ///
	+am_`y'[3,1]+am_`y'[3,2]+am_`y'[3,3]+am_`y'[3,4] ///
	+am_`y'[4,1]+am_`y'[4,2]+am_`y'[4,3]+am_`y'[4,4] ///
	+am_single_m_`y'[1,1]+am_single_m_`y'[2,1]+am_single_m_`y'[3,1]+am_single_m_`y'[4,1] ///
	+am_single_f_`y'[1,1]+am_single_f_`y'[1,2]+am_single_f_`y'[1,3]+am_single_f_`y'[1,4]
	
sca number_counterfactual=am_cf_sort_2018[1,1]+am_cf_sort_2018[1,2]+am_cf_sort_2018[1,3]+am_cf_sort_2018[1,4] ///
	+am_cf_sort_2018[2,1]+am_cf_sort_2018[2,2]+am_cf_sort_2018[2,3]+am_cf_sort_2018[2,4] ///
	+am_cf_sort_2018[3,1]+am_cf_sort_2018[3,2]+am_cf_sort_2018[3,3]+am_cf_sort_2018[3,4] ///
	+am_cf_sort_2018[4,1]+am_cf_sort_2018[4,2]+am_cf_sort_2018[4,3]+am_cf_sort_2018[4,4] ///
	+am_single_m_cf_sort_2018[1,1]+am_single_m_cf_sort_2018[2,1]+am_single_m_cf_sort_2018[3,1]+am_single_m_cf_sort_2018[4,1] ///
	+am_single_f_cf_sort_2018[1,1]+am_single_f_cf_sort_2018[1,2]+am_single_f_cf_sort_2018[1,3]+am_single_f_cf_sort_2018[1,4]
	
sca am_num_2=number_factual/(number_factual+number_counterfactual)

*denominator 2 /*The prob. of being in counterfactual sample*/
sca am_de_2=number_counterfactual/(number_factual+number_counterfactual)


*weights
mat am_w_couples=J(4,4,.)

forvalues i=1(1)4{
	forvalues j=1(1)4{
		mat am_w_couples[`i',`j']=(am_num_1_couples[`i',`j']/am_de_1_couples[`i',`j'])*(am_num_2/am_de_2)
	}
}

mat am_w_single_m=J(4,1,.)

forvalues i=1(1)4{
		mat am_w_single_m[`i',1]=(am_num_1_single_m[`i',1]/am_de_1_single_m[`i',1])*(am_num_2/am_de_2)
}

mat am_w_single_f=J(1,4,.)

forvalues i=1(1)4{
		mat am_w_single_f[1,`i']=(am_num_1_single_f[1,`i']/am_de_1_single_f[1,`i'])*(am_num_2/am_de_2)
}

gen weights=.
forvalues i=1(1)4{
	forvalues j=1(1)4{
replace weights=am_w_couples[`i',`j'] if am_male==`i' & am_female==`j' & relationship==1
	}
}

forvalues i=1(1)4{
	replace weights=am_w_single_m[`i',1] if relationship==0 & koen=="1" & ambition_type_k_4_s==`i'
}

forvalues i=1(1)4{
	replace weights=am_w_single_f[1,`i'] if relationship==0 & koen=="2" & ambition_type_k_4_s==`i'
}

replace weights=1 if weights==. /*Then those with missing type are counted same way in factual and counterfactual scenario*/

*Inequality 

qui ineqdeco household_earnings [weight=weights] if aar==`y'

sca gini_cf_sort_`y'=r(gini)
sca p90p50_cf_sort_`y'=r(p90p50)
sca p50p10_cf_sort_`y'=r(p50)/r(p10)



frame change results


replace gini_cf_sort=gini_cf_sort_`y' if aar==`y' 

replace p90p50_cf_sort=p90p50_cf_sort_`y' if aar==`y' 

replace p50p10_cf_sort=p50p10_cf_sort_`y' if aar==`y' 


frame change default


drop weights
}



*SAVE THE RESULTS*
frame change results

save "Results\tab_B3\Scenario_3\Decomposition\results_eq.dta", replace


